%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% O-Ring Production Networks
%   Author: Demir, Fieler, Xu, Yang
%   Last updated: January 2021
%
% Purpose: Standard Error (Full Model)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all; diary off; clc; clear; 


%% Parameters

    % load point estimates
    load('../../output/estimates/fullmodel_estimates.mat', 'est');            
    theta = est;

    % parameters
    define_param
    
    % increase upper bound
    param.ub = 10;
    
    
%% Variance-covariance matrix

    % load bootstrapped moments
    mmt_bootstrap_standard = xlsread('../../data/moments_bootstrap_standard.xlsx');
    mmt_bootstrap_clustered = xlsread('../../data/moments_bootstrap_clustered.xlsx');

    % variance-covariance matrix
    VC_standard = cov(mmt_bootstrap_standard);
    VC_clustered = cov(mmt_bootstrap_clustered);

    %check_diag = [diag(VC_standard)'; diag(VC_clustered)'];    
    
   
    
%% Jacobian (G)

    % initialize
    Nmom = size(VC_standard,1);
    Npar = size(theta,2);
    G = zeros(Nmom,Npar);
    iter_inner_h = zeros(1,Npar); iter_outer_h = zeros(1,Npar);
    iter_inner_l = zeros(1,Npar); iter_outer_l = zeros(1,Npar);
    
    % perturbation
    perturb = 0.05;
    h = abs(theta)*perturb;

    % jacobian
    for i = 1:Npar  
        
        fprintf('\nParameter %.0f\n', i);
        theta_h = theta; theta_h(i) = theta(i)+h(i);
        theta_l = theta; theta_l(i) = theta(i)-h(i);
        [moments_h, iter_h] = sim_moments(theta_h, param);
        [moments_l, iter_l] = sim_moments(theta_l, param);
        G(:,i) = (moments_h-moments_l)'./(2*h(i));        
        iter_inner_h(i) = iter_h(1); iter_outer_h(i) = iter_h(2);
        iter_inner_l(i) = iter_l(1); iter_outer_l(i) = iter_l(2);
        
    end    
    %save jacobian G iter_inner_h iter_outer_h iter_inner_l iter_outer_l
    %load jacobian   
    
    

%% Standard errors

    % number of simulations
    NS = param.Nf;

    % weighting matrix
    W_standard = inv(VC_standard);
    W_clustered = inv(VC_clustered);    
    
    % standard error
    Lambda_standard = (G'*W_standard*G)\(G'*W_standard);
    Lambda_clustered = (G'*W_clustered*G)\(G'*W_clustered);
    SE_standard = sqrt(diag((1+1/NS)*Lambda_standard*VC_standard*Lambda_standard'))';
    SE_clustered = sqrt(diag((1+1/NS)*Lambda_clustered*VC_clustered*Lambda_clustered'))';
    %save standard_error SE_standard SE_clustered

    % display
    %SE_standard
    %SE_clustered
    
        

%% Save Tables
diary off

    % create folder
    mkdir('../../output/tables') 

    % output 
    out_estimates_tables


